Frequency Tracking and Parameter Estimation for Robust Quantum State-Estimation 
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In this paper we consider the problem of tracking the state of a quantum system via a continuous 
measurement. If the system Hamiltonian is known precisely, this merely requires integrating the 
appropriate stochastic master equation. However, even a small error in the assumed Hamiltonian can 
render this approach useless. The natural answer to this problem is to include the parameters of the 
Hamiltonian as part of the estimation problem, and the full Bayesian solution to this task provides 
a state-estimate that is robust against uncertainties. However, this approach requires considerable 
computational overhead. Here we consider a single qubit in which the Hamiltonian contains a single 
unknown parameter. We show that classical frequency estimation techniques greatly reduce the 
computational overhead associated with Bayesian estimation and provide accurate estimates for the 
qubit frequency. 
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I. INTRODUCTION 



It is well known that the key to large scale quantum 
computing is the ability to maintain and manipulate a 
pure quantum state, minimizing any sources of possible 
error. Errors, either due to environmental fluctuations or 
inherent in the quantum system itself, are liable to pro- 
duce decoherence, dissipation or dephasing and lead to 
initially pure quantum states becoming mixed. In many 
of the candidate quantum technologies (such as quantum 
optics, atoms and ions in traps) the systems being used 
for as the qubits can be reasonably well characterised. 
Laser frequencies can be measured accurately. Atomic 
and ionic transition frequencies are consistent from one 
qubit to another. In the solid state, where devices need 
to be fabricated individually, there will always be some 
uncertainty in the parameters that determine the Hamil- 
tonians of the individual qubits, the coupling between 
qubits and the coupling to any external controls or envi- 
ronment. 

Any attempt to apply feedback control to solid-state 
qubits will therefore have to deal with uncertainties in 
their Hamiltonians. Feedback control is based upon con- 
tinuous state-estimation, in which the observer tracks the 
system via a continuous measurement P, [H, Q , and uses 
this knowledge to correct the evolution via control in- 
puts [1, H, H, 0| ■ It is therefore essential that continuous 
state-estimation handle uncertainties in the Hamiltonian. 
An estimation process that does this is referred to as 
being robust. An effective approach to achieving this ro- 
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bustness is to continuously estimate both the state of the 
system and the uncertain parameters in the Hamiltonian 
as the measurement proceeds. This is the approach we 
consider here - comparing the optimal Bayesian tech- 
nique and classical estimation techniques to determine 
an unknown parameter. 

Little work to date has examined continuous state- 
estimation in which the Hamiltonian is uncertain. Ex- 
ceptions are contained in work on parameter estimation 
by Gambetta and Wiseman 0, Verstraete et al. @ and 
Stockton et al. [To| . Gambetta and Wiseman examined 
a situation similar to the one we consider here, and de- 
rived the full Bayesian estimation equations. They found 
that ones ability to determine the unknown parameter 
depended strongly on way the system was measured. 
This is a result of the way in which the measured ob- 
servable is affected by the parameter, and the noise in- 
troduced into the dynamics by the measurement. The 
goal of the work by Verstraete et al. and Stockton et al. 
was to continuously estimate a signal that appeared as a 
(possibly time dependent) parameter in the Hamiltonian 
of a measured system. This problem arises naturally in 
the context of measuring classical fields, such as gravity- 
wave detection [ll[ and magnetometry [13]. In track- 
ing the classical parameter it was necessary to track the 
state of the system, and thus the authors also derived 
the full Bayesian estimation equations. More recently, 
Yamamoto [HI considered the problem of robust state- 
estimation for linear quantum systems, using results from 
classical control theory. Issues of robustness in quantum 
control have also been discussed in [TH, [l4[ . 

We consider what is possibly the simplest non-linear 
quantum state-estimation problem, that of tracking the 
state of a single qubit. We first examine what happens 
when the observer uses the usual stochastic master equa- 
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tion to estimate the state of the system, but with an 
incorrect Hamiltonian. We find that this procedure is, as 
expected, very sensitive to errors in the Hamiltonian, fail- 
ing quickly when the observer's Hamiltonian differs only 
a little from the true Hamiltonian. We then describe 
how this problem is solved by extending the estimation 
equations to include the Hamiltonian parameters. The 
resulting estimation equations are those derived in refer- 
ences [1] and [1]. (We also include a simple and concise 
derivation of these equations in the appendix.) Since 
the estimation equations simultaneously estimate both 
a quantum state and a set of classical parameters, we 
will refer to this set of equations as the hybrid stochastic 
master equation (HME). 

The HME, while an optimal solution to the robust 
state-estimation problem, does generate a considerable 
computational overhead. Our main purpose here is to 
show that this overhead can be reduced, with a rel- 
atively small reduction in performance, by using "fre- 
quency tracking" methods that were originally devel- 
oped for classical signal processing applications (such as 
radar frequency tracking, sonar processing and communi- 
cations) [ill. In particular, we examine three frequency 
tracking techniques: one Fourier-based method and two 
that do not require a Fourier transform (thereby avoid- 
ing a large computational overhead). The Fourier-based 
techniques is the maximiser of the periodogram - finding 
the frequency that maximises the correlation between the 
noisy measurement signal and a complex Fourier kcrnal. 
It is optimal for a single sinusoidal signal in noise, but 
it is computationally expensive. The first non-Fourier 
technique (Quinn-Fernandes [l6| ) is iterative and can be 
applied entirely in the time domain. It uses notch filter- 
ing and relies on the iterative construction of a filter that 
cancels an instability at the desired frequency. It can be 
shown to approach the Cramer-Rao bound for the conver- 
gence of the estimated frequency, in which the variance 
of the estimator asymptotes to 0(T~ 3 ), where T is the 
number of data points in the measurement record fl5j . 

The second non-Fourier technique is common in signal 
processing and is called 'multiple signal characterisation' 
or MUSIC [U [3 ■ This technique uses the properties of 
the eigenvectors of the autocovariance matrix of the data. 
These are related to the corresponding eigenvalues, which 
are in turn related to the frequencies present in the data. 
The variance of this technique scales as (9(T _1 ) rather 
than 0(T -3 ), but - unlike the Quinn-Fernandes method 
- it does not require an initial frequency estimate and can 
be more stable for large, appropriately filtered data sets. 
We compare the accuracy of these techniques with that 
of the full Bayesian estimator and find that they perform 
very well in general, although some do require additional 
filtering to obtain a stable parameter estimate. 

We begin by reviewing the formalism used to describe 
continuous weak measurement, and define the relevant 
measurement record upon which the parameter estimates 
will be based. We then describe the optimal Bayesian es- 
timator (the HME) for updating the parameter estimates 



and discuss several approximate methods for the estima- 
tion of frequencies from noisy data, with particular em- 
phasis on computational complexity. We then analyze 
the performance of the Fourier, Quinn-Fernandes and 
MUSIC techniques for a range of measurement strengths, 
and find that these perform well in comparison with the 
Bayesian method. We conclude by describing how these 
techniques can be included in a quantum feedback sys- 
tem to estimate an unknown (possibly time-dependent) 
Hamiltonian. 



II. CONTINUOUS WEAK MEASUREMENT 
AND STATE-ESTIMATION 

Measurement processes in quantum mechanics have 
always been the subject of intense theoretical interest. 
Traditionally, measurements have been projective and 
the quantum evolution has been reconstructed from a 
stochastic scries of results, averaged over an ensemble of 
similarly prepared quantum systems. The main problem 
with this approach is the restrictive nature of the projec- 
tive measurement process - there is no alternative to per- 
forming a series of projective measurements and building 
up a probabilistic model for the average evolution. Con- 
tinuous weak measurement models were developed as a 
means to describe the effect of weak (i.e. non-pro jective) 
measurement processes. This is achieved by coupling the 
quantum system of interest to environmental degrees of 
freedom, and then performing a series of projective mea- 
surements on the environment. Because the system be- 
comes correlated with the environment, the projective 
measurements provide a stream of information about the 
system. So long as the coupling strength between the 
system and the environment is sufficiently weak, then the 
effect of the projective measurements will be small com- 
pared to the underlying quantum evolution of the system, 
and the result is a (weak) continuous measurement. 

The evolution of the system now depends on the stream 
of measurement results, called the measurement record, 
and since these are necessarily random, the evolution of 
the system is given by a stochastic equation. If one aver- 
ages over all possible measurement records, then one ob- 
tains the usual Lindblad-form master equation for a sys- 
tem coupled to an environment. Each possible measure- 
ment record corresponds to a particular realisation of an 
experiment on a single quantum system. This is often re- 
ferred to as an 'unravelling' of the master equation, with 
different unravcllings corresponding to different kinds 
of measurements that one could make on the environ- 
ment [lj|. Weak measurements have a relatively small 
effect on the underlying quantum evolution, but increas- 
ing the strength of the coupling between the quantum 
system and the environment will also allow stronger mea- 
surements to be made - including projective-like mea- 
surements, where the non-unitary terms dominate the 
evolution - along with the quantum Zeno effect (20j . 

Of course, in a real experiment, all one has is the mea- 
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surcmcnt record rather than the underlying processes. 
In such cases it is possible to reconstruct aspects of 
the underlying quantum evolution by solving the rele- 
vant stochastic master equation (SME) for the specific 
measurement process/unravelling using the measurement 
record obtained from the experiment. The resultant evo- 
lution provides estimates of the instantaneous state of 
the quantum system given the measurement record [|J. 
This type of state estimation process is common in clas- 
sical systems, where the classical evolution of a system 
is estimated from a noisy series of measurements - often 
using the Kalman filter or its variants [2l[ . The parallels 
between classical state estimation and quantum state es- 
timation have been widely discussed in the literature [13] . 
In this context, the estimated quantum state represents 
the 'best guess' of what the actual state of the system is, 
based upon the set of measurements that are currently 
available. As the system evolves and more information 
becomes available via the measurements, the estimate of 
the quantum state should improve and the purity of the 
state should increase. 

For a quantum system with a Hamiltonian H, subject 
to a measurement corresponding to an operator c, the 
stochastic master equation for the unravelling is given 

by m, 



hdp c 



H, Pc 



1 



dt + I cp c tf - - (cc*p c + p c cc t ) J dt 



+dz^ (cp c - (c) Pc ) + (p c c t - Pc (cf)) dz (1) 



where dt is an infinitesimal time increment and dz is an 
infinitesimal complex Wiener increment. The complex 
Wiener increments obey the relations E[dz] = (where 
E[. . .} represents an expectation value), dzdz^ ~ HQdt, 
and dzdz T = KTdt. In general, dz is a complex vector 
and O and T are matrices, but if we restrict consider- 
ation to a single measurement interaction, is a real 
number (which we will set to one, corresponding to effi- 
cient detection (23|) and we set T = 1 so that the Wiener 
increments are real. p c is the density matrix conditioned 
upon the measurement record, which is given by [23[ 

dy{t) = (c T 6 + c t T) c dt + dz T = (c + c f ) c dt + dz (2) 

where the expectation values are taken with respect to 
the conditional state, (c) = Tr [p c c\. 

Unlike projective measurements, weak measurement 
processes are not restricted to Hermitian operators but 
- for the purposes of the current paper - we will assume 
that the measurement operator is Hermitian. This means 
that the sole action of the measurement is to extract in- 
formation about a single physical observable. (Measure- 
ment operators that are not Hermitian describe damping 
in addition to the extraction of information — a common 
example is the damping of a lossy optical cavity [HI)- For 
a Hermitian weak measurement c — \j2khy = , where k 
is the strength of the measurement interaction, the SME 



reduces to, 



dp c 



dt - k [y [y, p c }} dt 



+V2k(yp c + p c y-2(y) c p c )dW (3) 

where dW is a real Wiener increment (such that dW 2 = 
dt), and the measurement record is 



dy(t) 



= V8k(y) c dt + dW 



(4) 



At each point in time, the measurement record (from an 
experiment or a simulated process) will determine the as- 
sociated dW , given the current estimate of the state p c . 
The state of the system may then be updated using the 
value found for dW and equation ([5]). In classical state 
estimation, the value dW is often referred to as the inno- 
vation plj | , in that it represents the new information pro- 
vided by the measurement. The innovation is normally 
the difference between the actual measurement (in this 
case, dy(t)) and the current estimate of what the mea- 
surement should be given the current state of knowledge 
of the system (which is V&k {y) c dt). The main difference 
between classical and quantum measurement processes is 
that if the measured observable does not commute with 
the Hamiltonian, then the measurement actually feeds 
noise into the observable that is being measured. Be- 
cause of this, no matter how long we measure for, we will 
always be left with some uncertainty about the observ- 
able. 

The main potential problem with the above approach 
to state-estimation is that the stochastic master equation 
requires that the system is well-defined, in the sense that 
the Hamiltonian and the measurement interaction need 
to be known very accurately. If the Hamiltonian and in- 
teraction are not well defined, the measurement record 
that is used to update the system will produce incor- 
rect state estimates and assign too much confidence to 
the results of the SME. This is also common in classical 
state estimation problems, where too much confidence is 
assigned to the dynamical processes and valid measure- 
ments are often discounted because they conflict with the 
(inaccurate) state estimates [2l|. Ideally, we would like 
our state estimates to be insensitive to small errors in the 
Hamiltonian, or to allow the parameters in the Hamilto- 
nian to be estimated from the measurement record so 
that errors can be corrected on-line. 

To examine the sensitivity of the state estimates to 
errors in the Hamiltonian, we consider two SMEs: one 
is used to produce a measurement record (playing the 
role of an experimental system) and the other is used 
to estimate the underlying state of the system based on 
a Hamiltonian that contains a small error. Initially, we 
select a simple single-qubit Hamiltonian H = huj x a x /2 
(where o~ x is the Pauli x matrix), which generates a ro- 
tation of the single qubit Bloch vector about the Bloch 
x-axis. We then apply an interaction that generates a 
measurement along the Bloch z-axis. This maximises the 
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size of the output signal relative to the noise by choos- 
ing a measurement that is orthogonal to the axis of ro- 
tation. The Hamiltonian and measurement are selected 
because they reflect a common situation for solid state 
charge qubits, where the a x term corresponds to a tun- 
neling interaction and the a z interaction corresponds to 
a measurement of charge (see reference [25j for exam- 
ple) . Rewriting the SME in the Bloch vector representa- 
tion (Hi, with 

n = Tr[<7ip] i = X,Y, Z, 

where erj is a Pauli matix and the density operator is 

_ I + r x <J y + r Y <Ty + r z <j y 
P ~ 2 

the SME becomes three coupled stochastic equations 

drx = —4:krxdt — V8krxrzdW, 

dry = —uJ x Tzdt — Akrydt — V&kryrzdW, (5) 

dr z = u) x r Y dt + - r%)dW. 

The corresponding measurement record is, 

dy{t) 



8kr z dt + dW. 



(0) 



We use one SME with a fixed frequency u> x to generate 
a measurement record, and then feed this measurement 
record to a second SME. This second SME is used by an 
observer to estimate the state of the system. The first 
SME is initialised in a pure state. The second SME is ini- 
tialised in a completely mixed state [rx — r Y = Tz = 0), 
reflecting the fact that the observer does not initially 
know what state the qubit is in. If the Hamiltonian used 
in the second SME is perfect (that is, is exactly the same 
as that used in the first SME), then the initial mixed 
state will gradually purify as information is extracted. It 
will converge to the corrrect underlying state. (The pro- 
cess of purification - or state reduction - of single qubits 
subject to weak measurements and the use of feedback 
to increase the rate of purification have been discussed in 
detail by the current authors and others (2f| H3, HE ) • 
If the Hamiltonian used in the second SME contains an 
error, the estimated/conditioned state will purify, but 
the errors in the predicted evolution due to the inaccu- 
rate Hamiltonian will cause the estimated state and the 
underlying state from the first SME to diverge. The ques- 
tion is how accurate the Hamiltonian needs to be so that 
the system purifies before the prediction errors accumu- 
late and cause the states to diverge. A robust parame- 
ter estimation process would allow a system to be puri- 
fied to the required level before they diverge significantly. 
Preferably, this should be achieved without affecting the 
underlying evolution of the system, thereby allowing the 
Hamiltonian parameters to be estimated and updated. 
The alternative, a measurement strong enough to com- 
pletely dominate the evolution, projecting the state on a 



time-scale fast compared to the Hamiltonian, would cer- 
tainly give the correct estimated state, but it is not the 
situation we are interested in here. 

We begin by examining the divergence of the under- 
lying state, given by the first SME, from that of the es- 
timated state, given by the second SME. We do this by 
calculating the average fidelity between the two density 
matrices as a function of time for several different mea- 
surement strengths. The question relating to the extrac- 
tion of the parameters in the Hamiltonian will be dealt 
with in the next section. The fidelity F for two density 
matrices, po and p c , is a common measure of the differ- 
ence between two (possibly mixed) density matrices and 
it is given by [3C| . 



F = F(p Q ,p c ) = 



Tr 



In Fig. HH we show the average fidelity for a completely 
mixed initial state as a function of time for a num- 
ber of different expected frequency errors (the oscilla- 
tion frequency oj x being the parameter that we will be 
trying to estimate and the errors being characterised by 
their standard deviation as Ux - which will normally be 
given as a percentage) and for a measurement strength 
k = 0.07u x / (2ir). The fidelity is calculated by comparing 
the estimated density matrix p c and the underlying (pure 
state) density matrix pg. Both SMEs are integrated nu- 
merically using Milstein stochastic integration [3l| and at 
least 4000 time-steps per oscillation cycle. The estimated 
states purify as information is extracted from the mea- 
surement record, thereby increasing their fidelity. The 
average purity of p c will increase with time (the estimated 
purity tends to one) but the average fidelity eventually 
saturates as the information gain from the measurement 
record is balanced by the information loss due to the in- 
accuracies in the estimated Hamiltonian parameter. In 
the inset figure in Fig. 1, the saturation value of the aver- 
age fidelity is given as a function of measurement strength 
for three different values_of the error (as Ux = 1%, 2% and 
5%). Fig. 2 shows 1 — (F) (1— Average Fidelity, some- 
times called the 'infidelity') as a function of the percent- 
age error in the estimate of the oscillation frequency for 
three different measurement strengths. It is clear from 
the data shown in Figure 1 and Figure 2 that errors of 
the order of a few percent can have a significant effect 
on the fidelity of the estimation process. To obtain a fi- 
delity value in excess of 0.99 the frequency error must be 
less than about 0.2 — 0.3% if the measurement strength 
is around k = 0.03u; x /27T, but the error can be as large 
as 1% if the measurement strength is greater than about 
0.07 — 0.08oj. t /27t. Therefore, it is crucial that the mea- 
surement strength is chosen to reflect the expected ac- 
curacy of the frequency estimation process - although, 
the stronger the coupling between the measurement en- 
vironment and the qubit, the greater the effect of the 
back-action of the measurement on the qubit evolution. 
In the next section, we consider how an optimal Baycsian 
inference method can be used to update the parameter 
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FIG. 1: (Color online) Average fidelity of the estimated state 
vs. time for different values of expected errors asui^ = 1% 
(blue-solid), 2% (red-dash) and 5% (green-dot) for k = 
0.07^/(2-71-). Inset shows the saturated value for the aver- 
age fidelity as a function of measurement strength 2nk/u) x for 
P6u x = 1% (blue-solid), 2% (red-dash) and 5% (green-dot). 
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Fractional Error in Frequency 

FIG. 2: (Color online) Average fidelity error (1— Average Fi- 
delity) of the estimated state vs. percentage frequency er- 
rors for 2irk/u x = 0.03 (blue-solid), 0.07 (red-dash) and 0.12 
(green-dot) after 25 oscillation cycles. 



estimates on-line and consider several alternatives taken 
from classical signal processing and frequency tracking. 



III. QUANTUM PARAMETER ESTIMATION 
AND FREQUENCY TRACKING 

A. Bayesian Inference 

Making a continuous measurement will, in general, 
provide information not only about the state of a sys- 
tem, but also about the dynamics, and thus the parame- 



ters that determine the Hamiltonian. We will denote the 
parameters by the vector A. Our task is to process the 
measurement results in order to extract that information. 
The full solution to this problem is given by applying 
Bayesian inference [32| (classical measurement theory) 
in combination with quantum measurement theory (be- 
ing the quantum version of Bayesian inference [33|). To 
begin with, our state of knowledge regarding the system 
and the parameters is described by a density matrix for 
the state of the system, p c , and a probability density for 
the parameters, -P(A). Here we are concerned only with 
a single parameter (the frequency of oscillations) , so we 
have a scalar parameter, A. The solution to the continu- 
ous estimation problem tells us how to update both the 
density matrix and P(A) in each infinitesimal time-step, 
dt — that is, it gives differential equations for p c and 
P(A), and these are stochastic because they are driven 
by the measurement record. 

We now present the stochastic differential equations 
that solve the combined state and parameter estimation 
problem. While they look somewhat involved, they are 
quite simple to derive, and we give this derivation in Ap- 
pendix [XJ Since the Hamiltonian is a function of A, we 
will write it as H(X). The stochastic master equation for 
the density matrix is 



P(X)H(X)dX,p c 



dt - k[y, [y, p c ]]dt 



-Ak(yp c + p c y - 2(y) c p c )(dr - (y) c dt), (7) 



where dr is the measurement record, dr = dy(t)/y/8kh, 
y is the measured observable, and (y) c = Tr[yp c ]. The 
equation of motion for P(A) is 

dP(X) = 8k((y) - (y)x)(dr - (y)dt)P(X). (8) 

Here we have introduced a new quantity, (y)x- This is the 
expectation value that y would have, if we knew the value 
of A. To be able to calculate (y)x we must continually 
update our state-of-knowledge of the system, assuming 
that we know A, for every value of A. Denoting our state- 
of-knowledge of the system, given A, by p\, the update 
equation for each value of A is 



dp x = —[H(X),p x }}dt-k{y,{y,p x }}dt 

Px y~2(y) x p)(dr- (y) x dt) (9) 



h 

+4k{yp x 



where (y)\ = Tr[yp\]. We will refer to Eqs.flT])-© to- 
gether as the hybrid master equation, since it involves 
both quantum and classical states-of-knowlcdge. 

It is because we must propagate the density matrices 
px for every value of our parameter A that the full state- 
estimator is so numerically intensive. If A is a continuous 
variable, as is the case here, then it is impossible to have 
a density matrix for all values of A, and we therefore use 
a discrete grid of values. Note that the spacing of values 
on this grid must be a least as small as the accuracy with 
which we wish to track A. In addition, the length of the 
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grid (the range of values of A) must be large enough to 
encompass our initial uncertainty. If we use a fixed grid, 
these two requirements determine how many values of A 
require, which can be large. One approach to reducing 
this number would be to use a moving grid, updating the 
points on the grid along with -P(A). Even this procedure 
would require at least three copies of the density matrix 
(likely more), and certainly more if we had two or more 
parameters. In addition, it requires the rather involved 
task of interpolating a new density matrix each time a 
grid point is moved. Rather than exploring dynamical 
grids as a numerical method for implementing the full 
estimator, in the next section, we consider replacing the 
Bayesian updater for P(A) with a classical frequency es- 
timation technique. 

We present in Fig. [3] results for using the hybrid mas- 
ter equation to estimate the frequency of a single qubit. 
For these simulations we used a grid of 301 points for A. 
The simulations arc highly intensive: running the estima- 
tor for 500 periods, and averaging over 5000 realizations, 
takes ~ 33 hours on a 232 processor cluster. In Fig. Owe 
plot the average error of the estimator (the average rms 
deviation from the true frequency) as a function of time, 
and for different values of the measurement strength k. 
While a stronger measurement provides a faster rate of 
information extraction, the laws of quantum mechanics 
mean that it also introduces more noise into the system. 
As a result there is an optimal value for the measure- 
ment strength for a given frequency. We see from Fig. [3] 
that this is approximately k = 0.035w/27r, at least for a 
measurement time of 150 cycles. 
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FIG. 3: (Color online) Expected Frequency Errors (1 
standard deviation) vs. time for Bayesian Inference 
method and different measurement strengths, 2nk/uj x — 
0.005, 0.015, 0.025, 0.035, 0.040, 0.050. Inset figure shows the 
expected frequency error after 150 qubit cycles for the data 
shown in the main figure. 



B. Classical Frequency Tracking 

The ability to determine which frequencies are present 
in a noisy signal is extremely important in a large num- 
ber of practical engineering applications. In communi- 
cation systems, analogue signals are often transmitted 
by the modulation of the frequency of a carrier signal 
or digital data bits can be encoded in a frequency-shift 
key modulation scheme (34|. In classical signal process- 
ing, the frequencies present in sonar signals provide in- 
formation about the engine and gearing mechanisms in 
sea vessels [l5[ and tracking the Doppler frequency shifts 
in radar signals provides information about the relative 
(radial) motion of the radar and the target [35[ . A large 
number of techniques have been developed to extract the 
dominant frequencies from a noisy signal. Some rely on 
calculating the Fourier transform of the data record and 
estimating the likelihood that frequency component cor- 
responds to a signal rather than to noise, but these are 
computationally expensive and are similar in nature to 
the Bayesian method discussed above. However, other 
techniques have been developed which do not require ex- 
plicit calculation of a Fourier transform and some can be 
applied in the time domain alone. It is natural then to 
ask whether these techniques can be used to extract use- 
ful information from a quantum measurement record of 
the form given in ((6|) . From the wide variety of techniques 
that are available, we have chosen three for further study: 
the maximizer of the periodogram, the Quinn-Fcrnandes 
method [l6[ and the 'multiple signal characterisation' or 
MUSIC method [13, EH • The frequency that maximizes 
the periodogram provides the best correlation between 
the measurement record and a sine wave - it is a compu- 
tationally inefficient technique, but it provides the most 
accurate frequency estimate for cases where a single fre- 
quency signal is corrupted by white noise. The Quinn- 
Fcrnandes method is selected for its relative simplicity 
- it works in the time domain, removing the need to 
perform computationally expensive Fourier transforms - 
and for the fact that its theoretical accuracy approaches 
the optimal Cramer-Rao bound. The MUSIC method 
is slightly more complicated - using an estimate of the 
signal's autocorrelation matrix - but it does not require 
the initial frequency estimate used by Quinn-Fcrnandes 
and it can be used for signals containing more than one 
frequency component. 

The main advantage of classical techniques is that the 
measurement record can be integrated over much longer 
time steps, thereby reducing the computational require- 
ments. The optimal Bayesian inference technique re- 
quires that the SME is integrated using very small time 
steps. For the purity value of the quantum state to be 
well-behaved (i.e. not significantly larger than one), sev- 
eral thousand time steps are required per period of the os- 
cillation even when Milstein stochastic integration tech- 
niques are used. By contrast, the classical tracking tech- 
niques only require on the order of fifty sample points 
per period - a significant saving in computational terms. 
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The measurement record used by the classical trackers is 
a discrete sequence of values sampled at regular intervals 
from the continuous measurement record. As such, each 
sample, y n , is an integral of dy(t) over a finite time-step 
tn — > t„ + At where n = 0, ...TV — 1: 



*„+A* 



dy(t) 



(10) 



1. Maximizing the periodogram 

The periodogram for TV discrete points is formed by 
calculating 



JV-l 



y n exp(-itjt„) 



n=0 



(11) 



and the estimate of the qubit frequency is given by the 
value of omega that maximises I y {ui). For the cases con- 
sidered in this paper, the number of points will be large 
and the signal-to-noise ratio for the signal/measurement 
record will be very low when compared to most signal 
processing applications. Typically, there will be fifty 
points per oscillation cycle of the qubit, and we will con- 
sider the performance over several hundred cycles. The 
simple signal-to-noise ratio (simply the ratio between the 
two terms in the dy(t) measurement record) is given by: 



SNR = 10 log! 



(VSkAty 
2At 



101og 10 (4fcAt) 



which, for the discrete time examples that we will be 
considering, will provide signal-to-noise ratios between 
— 30dB and — 20dB. However, this simple calculation ig- 
nores the back-action of the measurement on the qubit, 
which induces slightly higher levels of noise in the output 
signal - and broadens the peak associated with the sinu- 
soidal oscillation of the qubit. In such cases, the accuracy 
of the Fouirer-bascd technique is limited by the broaden- 
ing of the peak because a broad peak has many maxima 
and a large number of points can generate a maximum 
when subject to a large amount of noise. Other tech- 
niques based on the Fourier transform - such as the in- 
terpolation of Fourier coefficients [36| - offer better per- 
formance than simply maximizing the discrete Fourier 
transform of the discrete time signal (which was found 
to provide very poor performance for all the cases stud- 
ied in this paper). However, these techniques will al- 
ways under-perform when compared to maximizing the 
periodogram. An alternative approach based on fitting 
the whole frequency response of the broadened measure- 
ment peak could be envisaged, and could offer more ac- 
curate frequency estimates, but this is likely to require 
even more computational power than calculating the pe- 
riodogram. 



2. The Quinn-Fernandes technique 

This technique relies on the fact that a sinusoidal sig- 
nal with additive noise will generate a frequency response 
which has a single peak and a noisy background [16[ . If 
a filter can be introduced to suppress the peak from the 
sinusoidal signal (a 'notch' filter) the remaining back- 
ground distribution will be the same as the underlying 
noise. This uses pole-zero cancellation. The filter is de- 
signed to have zero response at one frequency (the 'notch' 
frequency) and it is constructed to cancel the peak cor- 
responding to the sinusoidal frequency. If the signal is 
given by, 



A cos{nujAt + 



(12) 



for n — 0, ...N — 1, where A and <j> are the amplitude and 
phase of the sinusoidal signal and e„ is the noise, then 
the signal should satisfy the filter equation, 

x„-2cos(o;At)a;„_i+a;„_2 = e n -2cos(wAt)e n _i+e n _ 2 

(13) 

The Quinn-Fernandes algorithm provides a means to find 
the value of uj for which this equation is satisfied [l6l |. 
It operates purely in the time-domain and is therefore 
computationally efficient, over and above the efficiency 
savings associated with the low sampling rate used for 
the classical frequency estimation techniques used in this 
section. Rewriting Eq.(|13p as, 



(14) 



subject to the constraint a — /?, we start with an initial 
estimate u>\ for u, and set ot\ — 2 cos(wi At) along with an 
index j = 1. The data is then passed through a discrete- 
time filter to generate a new measurement record C, n j 



2 J 



(15) 



where ( n j = for n < 0. An estimate for j3 is then 
constructed, 



ft 



X)rt=0 (CnJ + Cn-2,j)Cn-l % 
Z-/n— ^n— l,j 



(16) 



If | ay — j3j | is larger than some threshold then the (3 value 
replaces the current a value, ay+i = /3j, and the pro- 
cess is applied again. If the difference between a and 
is small enough then the iteration in j is stopped and 
the signal frequency is given as u> = cos _1 (^/3j)/At. In 
most situations the number of iterations in j is fairly 
small, typically a few to ten iterations is sufficient to 
obtain a good frequency estimate, and the algorithm 
is also relatively insensitive to errors in the initial fre- 
quency estimate. However, these properties are nor- 
mally quoted for systems with higher signal-to-noise ra- 
tios than those considered here, and that do not gener- 
ate a back-action on the signal from the measurement 
process. In these slightly more benign situations the 



8 



Quinn-Fcrnandcs technique can be shown to generate er- 
rors which scale as 0(T~ 3 ) and approache the Cramer- 
Rao bound, as long as the error in the initial estimate is 
fairly small, OCT 1 - 1 ). 

For our discrete measurement record y n , we must pro- 
vide an initial estimate of the qubit frequency lo\ = Cu x . 
The accuracy of the initial estimate will have an effect 
on the accuracy of the final frequency estimate loq f ■ For 
sinusoidal signals with low to moderate noise levels the 
Quinn-Fernandcs technique is fairly insensitive to inaccu- 
racies in this initial estimate, but for the situation con- 
sidered here, where the signal-to-noise level is very low 
and the qubit oscillation is affected by the back-action 
of the measurement, the accuracy of the initial estimate 
is very important. Fig. 4 shows the evolution of the 
Quinn-Fernandes frequency estimates for the qubit sig- 
nal for an example measurement record, and for several 
different initial frequency estimates. On short timescales 
there are very few stable solutions to the iterative pro- 
cedure outlined above, but - as time increases - the ac- 
curacy of the initial estimate is required to be better be- 
cause the number of stable solutions/estimates increases. 
This multipicity of solutions is partly due to the very 
low signal-to-noise level of the measurement record and 
partly due to the broadening of the corresponding peak in 
the Fourier domain (one effect of measuring the qubit is 
to reduce the effective quality factor of the oscillation). 
The Quinn-Fernandes technique performs much better 
when the noise level is (artifically) lowered. Reducing the 
size of the dW term in the measurement record (equa- 
tion [5]) relative to the size of the signal reduces the num- 
ber of solutions, and the accuracy of the technique ap- 
proaches the theoretical best-case. (The best-case cor- 
responds to the output of the Quinn-Fernandes estimate 
when the initial estimate is perfect.) If the qubit sig- 
nal is a spin-i system, then moving to a larger value 
of angular momentum would boost the relative size of 
the signal compared to the dW noise and may give the 
Quinn-Fernandcs technique a significant advantage over 
the others. Increasing the signal-to- noise ratio by a factor 
of about 3 to 4 is enough to increase the accuracy of this 
technique significantly, and the Quinn-Fernandes tech- 
nique is the least expensive computationally of the dif- 
ferent techniques considered in this paper. Even without 
an increased signal-to-noise ratio, the Quinn-Fernandes 
does provide a slight improvement in accuracy over the 
initial estimate, as discussed in the next section. 



3. Multiple Signal Characterization (MUSIC) 

Multiple Signal Characterization, or MUSIC, is a 
widely used frequency estimation technique fUEl]- MU- 
SIC uses properties of the eigenvectors of the autocovari- 
ance matrix of the measurement record to decompose the 
data into a signal sub-space and a noise sub-space. The 
autocovariance matrix is formed by constructing a sym- 
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FIG. 4: (Color online) Example outputs from Quinn-Fernades 
frequency estimation algorithm as a function of time with 
slightly different initial frequency estimates. 
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(17) 



N-l 



c^n- 1 Y,{y n -y){y n -i-y) 



(18) 



where y is the mean of the measurement record. The 
size of the autocorrelation matrix M > 2 can be varied 
to optimize the performance of the technique and for the 
examples shown here M = 5, which provides a good level 
of performance and does not require large computational 
resources. MUSIC can be used to find more than one fre- 
quency, however, we will only use it to find one frequency 
in this paper. 

Once the autocovariance matrix has been constructed, 
the eigenvectors of the matrix v m (m = I...M) are found 
and ordered so that the eigenvalues 7 m are in ascending 
order 71 < 72 < . . . < 7m ■ If there are p sinusoids present 
in the data, the first M — p eigenvectors corresponding to 
the lowest M —p eigenvalues arc used to find the MUSIC 
'spectrum' S(co) from 



1 



Em=i |e*(wH 



(19) 



where e(w) = [1 exp(ia;Ai) cxp(2iwAt) ... cxp(i(M — 
1)ljA<)] t and the eigenvectors v m are normalized. The 
MUSIC spectrum contains a number of maxima, one for 
each of the p sinusoids present in the data. In our case, 
p = 1 and uj ~ oj x . The maximum value of the MUSIC 
spectrum therefore provides our estimate for the qubit 
frequency. 
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The signal-to-noise level in the measurement record is 
sufficiently high that the MUSIC technique applied on 
its own proves to be unstable - the frequency estimates 
generated occasionally contain very large errors, which 
unduly affect the variance of the errors and the resul- 
tant accuracy of the technique. To remove these very 
large errors, it is necessary to pre-filter the measurement 
record to remove frequency components too far away 
from the initial estimate of the qubit frequency. This 
is not as restrictive as the initial estimate used in the 
Quinn-Fcrnandes technique: the pre-filter used here is 
a bandpass filter (a Butterworth filter [37]]) with a pass 
band which is ±10% of the initial frequency estimate. 
Figure 5a shows an example of a MUSIC spectrum cal- 
culated for an example measurement record and the cor- 
responding MUSIC frequency estimates arc shown in Fig- 
ure 5b as a function of time. 




100 200 300 400 

Time (Number of Cycles) 

FIG. 5: (Color online) (a) Example MUSIC spectrum and (b) 
corresponding frequency estimate as a function of time. 



C. Comparison of Performance 

The principal aim of the frequency estimation algo- 
rithm is to improve the final fidelity of the estimated 



state in the stochastic master equation. Figures 1 and 
2 show that the expected fidelity of the state, after the 
initial purification phase, is dependent on the measure- 
ment strength and the expected error in the qubit fre- 
quency. The Bayesian method gives the best estimates 
for the qubit frequency for relatively weak measurements, 
2i:k/u! x ~ 0.035. However, such weak measurements re- 
quire very accurate frequency estimates to obtain high 
fidelity values: Figure 2 indicates that the frequency er- 
ror should be less than about 0.2 — 0.3% for the aver- 
age fidelity to be in excess of 0.99. Using stronger mea- 
surements, 2irk/u> x ~ 0.070 — 0.120, reduces the accu- 
racy required from the frequency estimation algorithm 
to around 0.8 — 1%. The strength of the measurement 
interaction is an important design consideration for solid 
state circuits because it is normally more difficult to de- 
couple or otherwise modify the interactions between fab- 
ricated circuits than it would be for an atomic or ionic 
qubit system. As such, the measurement strength would 
be an important design parameter in the construction of 
solid-state qubits and the associated control circuitry. 

The accuracy of the frequency estimation algorithms 
for the different measurement strengths is summarized 
in Figure 6. For weak measurement values around 
2irk/ui x ~ 0.030, the accuracy of each of the three classi- 
cal techniques is found to be approximately equal in the 
long time limit, although maximizing the periodogram 
is more accurate over short timescales. (Two values are 
given for the Quinn-Fcrnandes technique, representing 
the upper and lower bounds to the accuracy - which is 
a function of the accuracy of the initial estimate). The 
baseline accuracy of all three classical techniques is lim- 
ited to around 0.8%, which would be accurate enough 
for a maximum average fidelity of around F = 0.95 (see 
Figure 2). Only the Bayesian technique approaches the 
accuracy that would be required for fidelities approaching 
0.99 - at the expense of significantly increased computa- 
tional loading. 

As the measurement strength is increased, the accu- 
racy of the periodogram/Fourier-bascd technique is seen 
to deteriorate relative to either of the other two classical 
techniques. For very short periods of time, it remains the 
most accurate technique, but the accuracy asymptotes 
very quickly to an error of around 1.2%. With an initial 
1% error in the estimate, the Quinn-Fcrnandes technique 
provides a slight improvement in the frequency estimate 
after about 150 qubit cycles. However, the accuracy of 
the Quinn-Fernandes technique for 2irk/uj x ~ 0.070 is 
around 0.9%, which is just outside the desired range for 
an average fidelity of F = 0.99, and the accuracy of this 
approach can be seen to deteriorate when the measure- 
ment strength is increased. The lower bound for the 
Quinn-Fcrnandes accuracy is given in each graph in Fig- 
ure 6 and is the accuracy of the Quinn-Fernandes tech- 
nique if the initial estimate is perfect (i.e. zero error). 
This represents the best estimate that could be achieved 
using this (time-domain) technique, but in practice it 
will not be achievable because of errors in the initial frc- 
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qucncy. 

The expected accuracy of the MUSIC frequency esti- 
mates is slightly better than the other two classical meth- 
ods, the error reduces to around 0.7% at 2-Kkjuj x ~ 0.070. 
This is well within the range for an average fidelity 
greater than 0.99. The expected accuracy of this tech- 
nique would normally be expected to be worse than 
the other two methods because it does not converge 
to the Cramer- Rao bound in the long time limit, how- 
ever, it is found to be remarkably robust in the pres- 
ence of noise and it was also found to be stable in the 
presence of a broader peak in the Fourier domain. In 
fact, the expected accuracy was found to be fairly in- 
sensitive to noise and to variations in the measurement 
strength. It provides a reasonable estimate of the qubit 
frequency and is relatively efficient in terms of computa- 
tional demands - slightly longer to calculate than Quinn- 
Fernandes but much quicker that the pcriodogram or the 
Baycsian methods. 

Whilst the Baycsian frequency estimates are most ac- 
curate in the region of 2Ttk/u) x ~ 0.035, the accuracy of 
this technique after 500 cycles is always good enough to 
give an average fidelity above 0.99. However, comparing 
the trends for the fidelity errors shown in Figure 2 with 
the accuracies of the frequencies given in Figure 6, it sug- 
gests that the improvement in the average fidelity for in- 
creasing measurement strength increases faster than the 
degradation in the frequency estimates. This means that 
the Bayesian technique should still give higher fidelity 
estimates for the qubit states for stronger measurements 
even though frequency estimates are less accurate. It 
should also be noted that the Bayesian method generates 
a state estimate — to obtain an estimate of the state from 
any of the classical tracking methods one must integrate 
the SME once the qubit frequency has been determined. 




o.oi 
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IV. FREQUENCY TRACKING FOR A 
TIME-DEPENDENT HAMILTONIAN 

So far we have considered a frequency estimation pro- 
cess rather than a true tracking process, but the exten- 
sion of the estimation processes discussed to tracking is 
relatively straightforward for the case of a single qubit 
subject to a continuous measurement. The simplicity 
derives from the fact that there is only one frequency to 
be estimated in the measurement signal. If the qubit fre- 
quency is fixed, then more measurement data can be col- 
lected and the accuracy of the frequency estimates can be 
improved beyond the values given in Figure 6. The accu- 
racies of the Bayesian and MUSIC techniques are still im- 
proving after 500 qubit cycles. It is not clear whether this 
is also the case for the pcriodogram or Quinn-Fcrnandes 
techniques because any improvement in accuracy after a 
few hundred cycles is very slow in these cases. 

If the qubit frequency is not fixed but varies slowly 
with time - slow compared with the timescales associ- 
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FIG. 6: (Color online) Expected Frequency Errors (1 stan- 
dard deviation) vs. time for classical frequency tracking tech- 
niques: Bayesian (blue- (lower) solid), the maximum peri- 
odogram (red-cross), Lower limit for Quinn-Fernades (green- 
dot) error, Upper Limit for Quinn-Fernandes (green- (upper) 
solid) with 1% initial error (blue-dot), and MUSIC (purple- 
dash), for three different measurement strengths 2nk/uj x = 
0.03, 0.07 and 0.12 



ated with the frequency estimation process - there are 
a number of approaches that could be taken depending 
upon the nature of the variations in the frequency. For 
solid-state qubits, the most common problems are likely 
to be cither a gradual frequency drift or a discontinuous 
jump. A gradual frequency drift could arise because the 
qubit controls are not absolutely stable and give rise to 
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small drifts in the applied control fields leading to slight 
changes in the qubit frequency. Alternatively, a con- 
trol could contain infrequent, but sudden, small jumps 
in voltage which alter the bias point of the qubit circuit 
in a random way. 

For gradual frequency drifts, a natural method for gen- 
erating a time-dependent frequency 'track' would be to 
generate a moving window for the measurement data so 
that the frequency is always estimated over a fixed length 
window of data where the frequency is approximately 
stationary. The estimates are continually updated as the 
window is moved. The length of the window can be var- 
ied depending upon the anticipated time-dependence of 
the underlying signal. This can be computationally ex- 
pensive, the data having to be reprocessed each time the 
window moves, but it provides incremental updates to 
the frequency and generates a natural track history for 
the frequency. To reduce computational demands, the 
measurement data can be divided up into fixed blocks 
with finite length, but the frequency updates are then 
discrete measurements associated with each block and a 
separate tracking algorithm - such as a Kalman filter [2l[ 
- is required to generate a continuous track in-between 
the estimates/updates. For sudden discrete frequency 
changes, there are a number of techniques based on hid- 
den Markov models to generate the discontinuous tran- 
sitions between different frequency values [l5T ]. In each 
case, it is often helpful to have a model for the dynam- 
ics of the time-dependence of the system to reduce any 
residual errors in the tracked frequency values. 

The frequency estimation techniques can be general- 
ized if more than one frequency is present in the mea- 
surement signal. The MUSIC technique, in particular, is 
intended to allow multiple frequencies to be estimated. 
However, once measurements are made of more than 
one signal, potential problems can occur and slightly 
more thought is required when tracking multiple signals. 
When two frequencies are close together there are of- 
ten problems deciding which measurement is associated 
with which tracked frequency. This is a well-known prob- 
lem in object tracking. The mis-association of measure- 
ments with tracks can lead to significant errors in the 
parameter estimates. For algorithms like the Kalman fil- 
ter, additional information regarding the expected errors 
in the parameters being tracked (the frequencies in this 
case) propagates with the estimates. The expected er- 
rors can then be used to minimize the probability of mis- 
association by finding the most likely association of mea- 
surements with tracks based on a statistically-weighted 
distance (e.g. the Mahalanobis distance [H|). 



V. CONCLUSIONS 

We have considered a number of frequency estimation 
techniques that can be used with continuous quantum 
measurement processes. The techniques presented have 
been used to estimate the Hamiltonian for a single qubit, 



but they are quite general and can be readily adapted 
to monitor more complicated continuous measurement 
processes. The qubit has one parameter, the oscillation 
frequency for rotations around the Bloch x— axis. The 
most accurate estimates for the qubit frequency were gen- 
erated by a Bayesian inference technique, but this tech- 
nique is computationally expensive and typically requires 
thousands of integration steps per qubit cycle. By con- 
trast, the three classical frequency estimation techniques 
considered - maximized periodogram, Quinn-Fernandes, 
and MUSIC algorithms - all performed relatively well 
using only 50 integration steps per qubit cycle and one 
(MUSIC) was found to be robust to measurement noise 
and produced accurate frequency estimates for a wide 
range of measurement strengths. The classical frequency 
estimation algorithms are much less computationally in- 
tensive that Bayesian inference, partly due to the fact 
that they can be implemented with much larger time- 
steps. They provide frequency estimates that are ac- 
curate enough to generate robust (high fidelity) condi- 
tional state estimates for the qubit. We have also con- 
sidered how these frequency estimates can be integrated 
into a tracking system for a time-dependent Hamiltonian, 
where the qubit frequency undergoes continuous drifts or 
sudden discontinuous changes. 



APPENDIX A: DERIVING THE HYBRID 
MASTER EQUATION 

The observer starts with a state of knowledge defined 
by her density matrix for the system, p, and her proba- 
bility density for the parameters that define the Hamil- 
tonian of the system. For simplicity we will assume that 
there is only one real parameter A, so that the probability 
density is -P(A). Generalizing the result to any number 
of parameters is straightforward. 

If the observer knew A precisely, then upon obtain- 
ing the measurement result dr = (A)dt + dW/y/Bk, the 
evolution of the density matrix would be given by the 
stochastic master equation [l[ 

dpx = ^[H(X),p\dt-k[A,[A,p\}dt 

+4k(Ap + P A~2(A)p)(dr - (A)dt), (Al) 

so that the observers density matrix at time t + dt is 
px(t + dt) = p + dp\. Since the observer does not know 
which value A takes, she must average over all the possible 
p\ (t + dt) . Thus her state of knowledge at time t + dt is 
actually 

/oo 
P(X)p x (t + dt)dX 
-OO 

/OO 
[P(X)d Px ] dX 
-OO 

= p + dp. (A2) 



12 



The master equation for the observer's density matrix is 
therefore 



dp = T 



P{X)H(X)dX,p 



dt-k[A, [A,p]]dt 



+Ak{A P + P A- 2(A)p)(dr - (A)dt) (A3) 

We now have the equation that updates the observer's 
density matrix, but we need also the equation that up- 
dates the observers probability density P(A). To derive 
this we must use Bayes' theorem [38[. This tells us that 
the observer's new state of knowledge after obtaining the 
measurement result dr is 



P(X\dr) = P(dr\X)P(X)/Af, 



(A4) 



in which P(dr\X) is the probability density for dr given 
A, and Af is a normalization constant . If the value of the 
parameter was A, then the measurement record would be 
given by dr = (A) x dt + dW/V8k, where (A)\ = TrL4p A ], 
and p\ is obtained by using the update equation Eq. (|Al[) . 
Since dW is Gaussian distributed with mean and vari- 
ance dt, P(dr\X) is 



P(dr\X) = 



-4k(dr-(A) x dt) 2 /dt 



(A5) 



\/2ndt 

The updated probability density is thus 

P(X\dr) oc e -Mdr-(A) x dt) 2 /dtp^ ( A g) 

where dr is the actual measurement record obtained by 
the observer, dr = (A)dt + dW/VSk- We now expand 
this to first order in dt, remembering that dW 2 = dt [l[. 
This gives 

P(X\dr) OC e - ik ((A)dt-(A) x dt+dW/Vskf/dt p ^ 

e -4k((A) x -(A)) 2 dt+V8k((A) x -(Ay)dWp^ 

oc [l + V~8k((A) x - {A))dw\ P{X) (A7) 



The stochastic differential equation that updates the den- 
sity P(X) is therefore 



dP(X) 



V8k((A) x - (A))dWP(X) 

8k((A) x - (A))(dr - (A)dt)P(X). (A8) 



To be able to calculate this update to P(X) we must know 
(A}\ for all X. To track the quantum state along with the 
parameter A we must therefore solve the master equation 
for p (Eq. (|A2p ). the master equations for the p\ for every 
value of A (Eq. ljAljl ). and the update equation for P(X), 
Eq. (|A8|) , which is an example of a Kushner-Stratonovich 
equation. 
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